Statistical analysis

Litter decomposition

Source code and data found at https://github.com/maRce10/litter_decomposition_EFFEX

Purpose

• Evaluate role of nutrient availability on litter decomposition

 

Code

cols <- viridis(10, alpha = 0.7)
fill_color <- viridis(10)[7]

# brms models
chains <- 4
iters <- 40000
prior <- c(prior(normal(0, 10), "b"), prior(normal(0, 50), "Intercept"),
    prior(student_t(3, 0, 20), "sd"))

# set ggplot2 them
ggplot2::theme_set(theme_classic(base_size = 20))


# standard error
se <- function(x) sd(x)/sqrt(length(x))

1 Nutrient change

1.1 Nitrogen

Content

Code

nutr <- read.csv("./data/raw/litter-nutrients-mas.csv")

nutr$plot.f <- as.factor(nutr$plot)

nutr$days.sc <- scale(nutr$days)

nutr$litter.n.content.prop.initial <- nutr$litter.n.content.perc.initial/100

# remove plot 9
sub.nutr <- nutr[nutr$plot != 9, ]

agg_n <- aggregate(litter.n.content.perc.initial ~ colecta + treat +
    days, nutr, mean)

agg_n$sd <- aggregate(litter.n.content.perc.initial ~ colecta + treat +
    days, nutr, sd)$litter.n.content.perc.initial

agg_n$se <- aggregate(litter.n.content.perc.initial ~ colecta + treat +
    days, nutr, se)$litter.n.content.perc.initial

agg_n$treat <- factor(agg_n$treat, levels = c("C", "N", "P", "NP"))

pd <- position_dodge(15)

ggplot(agg_n, aes(x = days, y = litter.n.content.perc.initial, color = treat)) +
    geom_point(size = 2, position = pd) + geom_errorbar(aes(ymax = litter.n.content.perc.initial +
    se, ymin = litter.n.content.perc.initial - se), width = 0, position = pd) +
    geom_line(size = 1.2, position = pd) + scale_color_viridis_d(alpha = 0.5) +
    labs(x = "Time (days)", y = "Litter N content (% initial)", color = "Treatment") +
    scale_x_continuous(breaks = unique(agg_n$days), labels = unique(agg_n$days)) +
    theme(legend.position = c(0.9, 0.8))

Code

ggsave("./output/Litter_N_content.tiff", dpi = 300)

Download image

Code

fit.n <- brm(litter.n.content.prop.initial ~ treat * days.sc + (1 |
    plot.f), data = nutr, chains = chains, cores = chains, family = Beta(),
    iter = iters, control = list(adapt_delta = 0.99, max_treedepth = 15),
    prior = prior, file = "./data/processed/litter.n.content.rem_model",
    file_refit = "on_change")

Code

extended_summary(read.file = "./data/processed/litter.n.content.rem_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

1.2 litter.n.content.rem_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20)	litter.n.content.prop.initial ~ treat * days.sc + (1 | plot.f)	20000	4	1	10000	0	0	19631.06	22594.55	272268972

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
treatN	-0.340	-0.634	-0.048	1	27796.92	26514.82
treatNP	-0.051	-0.344	0.238	1	27673.41	24838.46
treatP	-0.060	-0.349	0.231	1	28095.99	25183.40
days.sc	-0.503	-0.712	-0.307	1	19631.06	22594.55
treatN:days.sc	-0.094	-0.369	0.183	1	25338.01	27759.09
treatNP:days.sc	0.039	-0.227	0.314	1	25217.60	27081.45
treatP:days.sc	0.048	-0.220	0.323	1	23691.70	28287.99

Concentration

Code

agg_conc_n <- aggregate(perc.n ~ colecta + treat + days, nutr, mean)

agg_conc_n$sd <- aggregate(perc.n ~ colecta + treat + days, nutr,
    sd)$perc.n

agg_conc_n$se <- aggregate(perc.n ~ colecta + treat + days, nutr,
    se)$perc.n

agg_conc_n$treat <- factor(agg_conc_n$treat, levels = c("C", "N",
    "P", "NP"))

pd <- position_dodge(15)

ggplot(agg_conc_n, aes(x = days, y = perc.n, color = treat)) + geom_point(size = 2,
    position = pd) + geom_errorbar(aes(ymax = perc.n + se, ymin = perc.n -
    se), width = 0, position = pd) + geom_line(size = 1.2, position = pd) +
    scale_color_viridis_d(alpha = 0.5) + labs(x = "Time (days)", y = "Litter N concentration (%)",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_conc_n$days),
    labels = unique(agg_conc_n$days)) + theme(legend.position = c(0.3,
    0.8))

Code

ggsave("./output/Litter_N_concentration.tiff", dpi = 300)

Download image

Code

nutr$prop.n <- nutr$perc.n/100

fit.perc.n <- brm(prop.n ~ treat * days.sc + (1 | plot.f), data = nutr,
    chains = chains, core = chains, family = Beta(), iter = iters,
    control = list(adapt_delta = 0.99, max_treedepth = 15), prior = prior,
    file = "./data/processed/litter.n.perc_model", file_refit = "on_change")

Code

extended_summary(read.file = "./data/processed/litter.n.perc_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

1.3 litter.n.perc_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20)	prop.n ~ treat * days.sc + (1 | plot.f)	20000	4	1	10000	0	0	16686.78	23193.79	757288462

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
treatN	-0.097	-0.240	0.046	1	25471.79	25380.38
treatNP	-0.052	-0.192	0.087	1	25656.85	26419.36
treatP	-0.025	-0.163	0.114	1	24461.38	24258.37
days.sc	0.134	0.041	0.227	1	16686.78	23193.79
treatN:days.sc	-0.075	-0.205	0.052	1	20777.70	25856.35
treatNP:days.sc	-0.030	-0.160	0.099	1	20665.11	26156.06
treatP:days.sc	-0.012	-0.138	0.116	1	20899.42	26866.06

1.4 Phosphorus

Content

Code

nutr$litter.p.content.prop.initial <- nutr$litter.p.content.perc.initial/100

# excluding plot 9
agg_p <- aggregate(litter.p.content.perc.initial ~ colecta + treat +
    days, sub.nutr, mean)

agg_p$sd <- aggregate(litter.p.content.perc.initial ~ colecta + treat +
    days, sub.nutr, sd)$litter.p.content.perc.initial

agg_p$se <- aggregate(litter.p.content.perc.initial ~ colecta + treat +
    days, sub.nutr, se)$litter.p.content.perc.initial


agg_p$treat <- factor(agg_p$treat, levels = c("C", "N", "P", "NP"))


pd <- position_dodge(15)

ggplot(agg_p, aes(x = days, y = litter.p.content.perc.initial, color = treat)) +
    geom_point(size = 2, position = pd) + geom_errorbar(aes(ymax = litter.p.content.perc.initial +
    se, ymin = litter.p.content.perc.initial - se), width = 0, position = pd) +
    geom_line(size = 1.2, position = pd) + scale_color_viridis_d(alpha = 0.5) +
    labs(x = "Time (days)", y = "Litter P content (% initial)", color = "Treatment") +
    scale_x_continuous(breaks = unique(agg_p$days), labels = unique(agg_p$days)) +
    theme(legend.position = c(0.9, 0.8))

Code

ggsave("./output/Litter_P_content.tiff", dpi = 300)

Download image

Code

fit.p <- brm(litter.p.content.prop.initial ~ treat * days.sc + (1 |
    plot.f), data = nutr, chains = chains, prior = prior, cores = chains,
    iter = iters, control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "./data/processed/litter.p.content.rem_model", file_refit = "on_change")

Code

extended_summary(read.file = "./data/processed/litter.p.content.rem_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

1.5 litter.p.content.rem_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	litter.p.content.prop.initial ~ treat * days.sc + (1 | plot.f)	20000	4	1	10000	0	0	13618.61	17352.52	2063397055

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
treatN	-0.090	-1.258	1.081	1	13662.76	18828.78
treatNP	0.009	-1.148	1.172	1	13618.61	17352.52
treatP	0.754	-0.399	1.907	1	14549.30	18770.67
days.sc	-0.134	-0.473	0.205	1	17312.08	23553.03
treatN:days.sc	0.001	-0.459	0.459	1	20983.83	27133.63
treatNP:days.sc	0.021	-0.441	0.482	1	21015.75	25715.30
treatP:days.sc	-0.168	-0.620	0.288	1	20970.54	26482.58

Concentration

Code

agg_conc_p <- aggregate(perc.p ~ colecta + treat + days, sub.nutr,
    mean)

agg_conc_p$sd <- aggregate(perc.p ~ colecta + treat + days, sub.nutr,
    sd)$perc.p

agg_conc_p$se <- aggregate(perc.p ~ colecta + treat + days, sub.nutr,
    se)$perc.p

agg_conc_p$treat <- factor(agg_conc_p$treat, levels = c("C", "N",
    "P", "NP"))

pd <- position_dodge(15)

ggplot(agg_conc_p, aes(x = days, y = perc.p, color = treat)) + geom_point(size = 2,
    position = pd) + geom_errorbar(aes(ymax = perc.p + se, ymin = perc.p -
    se), width = 0, position = pd) + geom_line(size = 1.2, position = pd) +
    scale_color_viridis_d(alpha = 0.5) + labs(x = "Time (days)", y = "Litter P concentration (%)",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_conc_p$days),
    labels = unique(agg_conc_p$days)) + theme(legend.position = c(0.9,
    0.8))

Code

ggsave("./output/Litter_P_concentration.tiff", dpi = 300)

Download image

Code

# convert to proportions to use beta distribution
sub.nutr$prop.p <- sub.nutr$perc.p/100

fit.perc.p <- brm(prop.p ~ treat * days.sc + (1 | plot.f), data = sub.nutr,
    chains = chains, cores = chains, family = Beta(), iter = iters,
    control = list(adapt_delta = 0.99, max_treedepth = 15), prior = prior,
    file = "./data/processed/litter.p.perc_model", file_refit = "on_change")

Code

extended_summary(read.file = "./data/processed/litter.p.perc_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

1.6 litter.p.perc_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20)	prop.p ~ treat * days.sc + (1 | plot.f)	20000	4	1	10000	0	0	17993.25	23361.17	1797840914

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
treatN	-0.106	-0.634	0.409	1	28569.54	28049.17
treatNP	0.031	-0.475	0.530	1	27422.98	28348.78
treatP	0.077	-0.468	0.610	1	28087.50	28163.62
days.sc	-0.034	-0.388	0.302	1	17993.25	23361.17
treatN:days.sc	-0.001	-0.486	0.478	1	22840.20	28837.19
treatNP:days.sc	0.040	-0.430	0.503	1	21964.87	27252.32
treatP:days.sc	0.028	-0.475	0.531	1	23659.98	28583.96

Takeaways

• Litter P content is significantly lower in plus N treatment plots than in control plot, after accounting for variation explained by time

 

1.7 Remaining litter

Download image

Code

fit <- brm(prop.litter.rem ~ trat * days.sc + (1 | plot.f), data = dat,
    chains = chains, family = Beta(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/prop.litter.rem_model",
    file_refit = "on_change")

# dat$days.fc <- as.numeric(as.factor(dat$days)) # monotonic
# effect of time fit_mo <- brm(prop.litter.rem ~ trat *
# mo(days.fc) + (1 | plot.f), data = dat, chains = chains,
# family = Beta(), iter = iters, control =
# list(adapt_delta=0.99, max_treedepth=15), cores = chains,
# prior = prior, file =
# './data/processed/prop.litter.rem_model_monotonic', file_refit
# = 'on_change')

Code

extended_summary(read.file = "./data/processed/prop.litter.rem_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

1.8 prop.litter.rem_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20)	prop.litter.rem ~ trat * days.sc + (1 | plot.f)	20000	4	1	10000	0	0	17979.85	23910.9	1631599382

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
tratN	-0.038	-0.203	0.125	1	30540.12	29362.46
tratNP	0.074	-0.085	0.233	1	29883.05	29157.97
tratP	0.032	-0.127	0.192	1	30065.88	30156.88
days.sc	-0.683	-0.800	-0.570	1	17979.85	23910.90
tratN:days.sc	-0.047	-0.207	0.113	1	23277.78	30456.01
tratNP:days.sc	0.106	-0.045	0.259	1	22364.84	28848.92
tratP:days.sc	0.056	-0.097	0.212	1	22771.51	27130.32

1.9 Remaining wood

Download image

 

Code

dat$prop.wood.rem <- dat$perc.wood.rem/100

fit2 <- brm(prop.wood.rem ~ trat * days.sc + (1 | plot.f), data = dat,
    chains = chains, family = Beta(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/prop.wood.rem_model",
    file_refit = "on_change")

Code

extended_summary(read.file = "./data/processed/prop.wood.rem_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

1.10 prop.wood.rem_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20)	prop.wood.rem ~ trat * days.sc + (1 | plot.f)	20000	4	1	10000	0	0	16872.71	23764.8	1828689869

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
tratN	0.118	-0.128	0.366	1	25420.93	24370.70
tratNP	0.254	0.012	0.496	1	24798.43	26473.85
tratP	0.197	-0.045	0.438	1	25725.64	25244.62
days.sc	-0.910	-1.094	-0.737	1	16872.71	23764.80
tratN:days.sc	-0.006	-0.244	0.234	1	20661.18	27810.64
tratNP:days.sc	0.097	-0.140	0.333	1	20685.17	28019.73
tratP:days.sc	0.175	-0.058	0.412	1	20811.62	27760.94

2 K

2.1 Litter

Correlation between Silvia’s and Andrea’s K

Code

k_vals <- read.csv("./data/raw/k-values-corr.csv")

cor(k_vals$av.k.litter, k_vals$sil.k.litter)

[1] 0.9327804

Code

cor(k_vals$av.k.wood, k_vals$sil.k.wood)

[1] 0.9621473

Comparing all treatments vs control

Code

agg_kvals <- aggregate(sil.k.litter ~ Treatment, k_vals, mean)
agg_kvals$se <- aggregate(sil.k.litter ~ Treatment, k_vals, se)$sil.k.litter
agg_kvals$sd <- aggregate(sil.k.litter ~ Treatment, k_vals, sd)$sil.k.litter
agg_kvals$n <- aggregate(sil.k.litter ~ Treatment, k_vals, length)$sil.k.litter
agg_kvals$treat <- factor(agg_kvals$Treatment, levels = c("C", "N", "P", "NP"))
agg_kvals$n.labels <- paste("n =", agg_kvals$n) 

# composed box plot
ggplot(k_vals, aes(x = Treatment, y = sil.k.litter)) +
## add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    ## custom bandwidth
    adjust = .5,
    ## adjust height
    width = .6,
    .width = 0,
    ## move geom to the right
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    ## remove outliers
    outlier.shape = NA ## `outlier.shape = NA` works as well
  ) +
  ## add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    ## draw jitter on the left
    side = "l",
    ## control range of jitter
    range_scale = .4,
    ## add some transparency
    alpha = .5,
    transformation = ggplot2::position_jitter(height = 0)

  ) +
   labs(y = "Decomposition constant (k)") +
  # ylim(c(-0.39, 0.145)) +
  geom_text(data = agg_kvals, aes(y = rep(0.4, nrow(agg_kvals)), x = Treatment, label = n.labels), nudge_x = 0, size = 6) +
     theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
    scale_x_discrete(labels=c("C" = "Control", "N" = "+N",
                              "NP" = "+NP", "P" = "+P"))

Code

ggsave("./output/litter_decomposition_k_by_treatment.tiff", dpi = 300)

Download image

Code

fit_k.litter <- brm(sil.k.litter ~ Treatment + (1 | quadrat), data = k_vals,
    chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/k_litter_model",
    file_refit = "on_change")

Code

extended_summary(gsub.pattern = "b_treatment|b_", gsub.replacement = "",
    highlight = TRUE, remove.intercepts = TRUE, read.file = "./data/processed/k_litter_model.rds")

2.2 k_litter_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.litter ~ Treatment + (1 | quadrat)	20000	4	1	10000	5	0	22486.82	24586.15	224004760

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
TreatmentN	0.238	-0.176	0.659	1	22486.82	24586.15
TreatmentNP	-0.156	-0.572	0.261	1	25034.00	26992.34
TreatmentP	-0.006	-0.425	0.412	1	23833.14	25611.23

Phosphorus vs no-phosphorus

Code

fit_k.litter.p.np <- brm(sil.k.litter ~ p.treat + (1 | quadrat), data = k_vals,
    chains = chains, family = gaussian(), iter = iters, , control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/k_litter_p_nop_model",
    file_refit = "on_change")

Code

extended_summary(read.file = "./data/processed/k_litter_p_nop_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

2.3 k_litter_p_nop_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.litter ~ p.treat + (1 | quadrat)	20000	4	1	10000	8	0	38432.94	27572.89	1236972411

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
p.treatwith.p	-0.201	-0.492	0.092	1	38432.94	27572.89

Code

agg_kvals <- aggregate(sil.k.litter ~ p.treat, k_vals, mean)
agg_kvals$se <- aggregate(sil.k.litter ~ p.treat, k_vals, se)$sil.k.litter
agg_kvals$sd <- aggregate(sil.k.litter ~ p.treat, k_vals, sd)$sil.k.litter
agg_kvals$p.treat <- factor(agg_kvals$p.treat, labels = c("No P", "P"))
agg_kvals$n <- aggregate(sil.k.litter ~ p.treat, k_vals, length)$sil.k.litter
agg_kvals$n.labels <- paste("n =", agg_kvals$n) 
k_vals$p.treat <- factor(k_vals$p.treat, labels = c("No P", "P"))

# composed box plot
rn_p_litter <- ggplot(k_vals, aes(x = p.treat, y = sil.k.litter)) +
## add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    ## custom bandwidth
    adjust = .5,
    ## adjust height
    width = .6,
    .width = 0,
    ## move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    ## remove outliers
    outlier.shape = NA ## `outlier.shape = NA` works as well
  ) +
  ## add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    ## draw jitter on the left
    side = "l",
    ## control range of jitter
    range_scale = .4,
    ## add some transparency
    alpha = .5,
    transformation = ggplot2::position_jitter(height = 0)

  ) +
  labs(y = "Decomposition constant (k)", x = "P treatment") +
  geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = p.treat, label = n.labels), nudge_x = 0, size = 6) +
     theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1))  +
    scale_x_discrete(labels=c("No P" = "No P", "P" = "+P"))

rn_p_litter

Code

ggsave("./output/litter_phosphorus_decomposition_k.tiff", dpi = 300)

Nitrogen vs no-nitrogen

Code

fit_k.litter.n.nn <- brm(sil.k.litter ~ n.treat + (1 | quadrat), data = k_vals,
    chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/k_litter_n_no_n_model",
    file_refit = "on_change")

Code

extended_summary(read.file = "./data/processed/k_litter_n_no_n_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

2.4 k_litter_n_no_n_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.litter ~ n.treat + (1 | quadrat)	20000	4	1	10000	17	0	37023.21	27505.83	956380015

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
n.treatwith.n	0.045	-0.254	0.345	1	37023.21	27505.83

Code

agg_kvals <- aggregate(sil.k.litter ~ n.treat, k_vals, mean)
agg_kvals$se <- aggregate(sil.k.litter ~ n.treat, k_vals, se)$sil.k.litter
agg_kvals$sd <- aggregate(sil.k.litter ~ n.treat, k_vals, sd)$sil.k.litter
agg_kvals$n.treat <- factor(agg_kvals$n.treat, labels = c("No N", "N"))
agg_kvals$n <- aggregate(sil.k.litter ~ n.treat, k_vals, length)$sil.k.litter
agg_kvals$n.labels <- paste("n =", agg_kvals$n) 
k_vals$n.treat <- factor(k_vals$n.treat, labels = c("No N", "N"))

# composed box plot
rn_n_litter <- ggplot(k_vals, aes(x = n.treat, y = sil.k.litter)) +
## add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    ## custom bandwidth
    adjust = .5,
    ## adjust height
    width = .6,
    .width = 0,
    ## move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    ## remove outliers
    outlier.shape = NA ## `outlier.shape = NA` works as well
  ) +
  ## add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    ## draw jitter on the left
    side = "l",
    ## control range of jitter
    range_scale = .4,
    ## add some transparency
    alpha = .5,
    transformation = ggplot2::position_jitter(height = 0)

  ) +
  labs(y = "Decomposition constant (k)", x = "P treatment") +
  geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = n.treat, label = n.labels), nudge_x = 0, size = 6) +
     theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1))  +
    scale_x_discrete(labels=c("No P" = "No P", "P" = "+P"))

rn_n_litter

Code

ggsave("./output/litter_nitrogen_decomposition_k.tiff", dpi = 300)

2.5 Wood

Comparing all treatments vs control

Code

agg_kvals <- aggregate(sil.k.wood ~ Treatment, k_vals, mean)
agg_kvals$se <- aggregate(sil.k.wood ~ Treatment, k_vals, se)$sil.k.wood
agg_kvals$sd <- aggregate(sil.k.wood ~ Treatment, k_vals, sd)$sil.k.wood
agg_kvals$n <- aggregate(sil.k.wood ~ Treatment, k_vals, length)$sil.k.wood
agg_kvals$treat <- factor(agg_kvals$Treatment, levels = c("C", "N", "P", "NP"))
agg_kvals$n.labels <- paste("n =", agg_kvals$n) 

# composed box plot
ggplot(k_vals, aes(x = Treatment, y = sil.k.wood)) +
## add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    ## custom bandwidth
    adjust = .5,
    ## adjust height
    width = .6,
    .width = 0,
    ## move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    ## remove outliers
    outlier.shape = NA ## `outlier.shape = NA` works as well
  ) +
  ## add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    ## draw jitter on the left
    side = "l",
    ## control range of jitter
    range_scale = .4,
    ## add some transparency
    alpha = .5,
    transformation = ggplot2::position_jitter(height = 0)

  ) +
   labs(y = "Decomposition constant (k)") +
  # ylim(c(-0.39, 0.145)) +
  geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = Treatment, label = n.labels), nudge_x = 0, size = 6) +
     theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
    scale_x_discrete(labels=c("C" = "Control", "N" = "+N",
                              "NP" = "+NP", "P" = "+P"))

Code

ggsave("./output/wood_decomposition_k_by_treatment.tiff", dpi = 300)

Download image

Code

fit_k.wood.p.treat <- brm(sil.k.wood ~ Treatment + (1 | quadrat),
    data = k_vals, chains = chains, family = gaussian(), iter = iters *
        2, control = list(adapt_delta = 0.99, max_treedepth = 15),
    cores = chains, prior = prior, file = "./data/processed/k_wood_treatment_model",
    file_refit = "on_change")

Code

extended_summary(read.file = "./data/processed/k_wood_treatment_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

2.6 k_wood_treatment_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.wood ~ Treatment + (1 | quadrat)	80000	4	1	40000	23	0	96744.72	104808.8	528819192

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
TreatmentN	0.034	-0.537	0.610	1	96900.55	107769.7
TreatmentNP	-0.490	-1.064	0.087	1	96744.72	104808.8
TreatmentP	-0.469	-1.044	0.105	1	97055.26	107363.6

3 Phosphorus vs no-phosphorus

Code

agg_kvals <- aggregate(sil.k.wood ~ p.treat, k_vals, mean)
agg_kvals$se <- aggregate(sil.k.wood ~ p.treat, k_vals, se)$sil.k.wood
agg_kvals$sd <- aggregate(sil.k.wood ~ p.treat, k_vals, sd)$sil.k.wood


agg_kvals$p.treat <- factor(agg_kvals$p.treat, labels = c("No P", "P"))
agg_kvals$n <- aggregate(sil.k.wood ~ p.treat, k_vals, length)$sil.k.wood
agg_kvals$n.labels <- paste("n =", agg_kvals$n) 


k_vals$p.treat <- factor(k_vals$p.treat, labels = c("No P", "P"))


# composed box plot
rn_p_wood <- ggplot(k_vals, aes(x = p.treat, y = sil.k.wood)) +
## add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    ## custom bandwidth
    adjust = .5,
    ## adjust height
    width = .6,
    .width = 0,
    ## move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    ## remove outliers
    outlier.shape = NA ## `outlier.shape = NA` works as well
  ) +
  ## add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    ## draw jitter on the left
    side = "l",
    ## control range of jitter
    range_scale = .4,
    ## add some transparency
    alpha = .5,
    transformation = ggplot2::position_jitter(height = 0)

  ) +
  labs(y = "Decomposition constant (k)", x = "P treatment") +
  geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = p.treat, label = n.labels), nudge_x = 0, size = 6) +
     theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1))  +
    scale_x_discrete(labels=c("No P" = "No P", "P" = "+P"))

rn_p_wood

Code

ggsave("./output/phosphorus_decomposition_k.tiff", dpi = 300)

Download image

Code

fit_k.wood.p.no.p <- brm(sil.k.wood ~ p.treat + (1 | quadrat), data = k_vals,
    chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/k_wood_p_no_p_model",
    file_refit = "on_change")

Code

extended_summary(read.file = "./data/processed/k_wood_p_no_p_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

3.1 k_wood_p_no_p_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.wood ~ p.treat + (1 | quadrat)	20000	4	1	10000	4	0	31952.25	26238.54	326631910

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
p.treatwith.p	-0.496	-0.892	-0.1	1	31952.25	26238.54

4 Nitrogen vs no-nitrogen

Code

agg_kvals <- aggregate(sil.k.wood ~ n.treat, k_vals, mean)
agg_kvals$se <- aggregate(sil.k.wood ~ n.treat, k_vals, se)$sil.k.wood
agg_kvals$sd <- aggregate(sil.k.wood ~ n.treat, k_vals, sd)$sil.k.wood
agg_kvals$n.treat <- factor(agg_kvals$n.treat, labels = c("No N", "N"))
agg_kvals$n <- aggregate(sil.k.wood ~ n.treat, k_vals, length)$sil.k.wood
agg_kvals$n.labels <- paste("n =", agg_kvals$n) 
k_vals$n.treat <- factor(k_vals$n.treat, labels = c("No N", "N"))

# composed box plot
rn_n_wood <- ggplot(k_vals, aes(x = n.treat, y = sil.k.wood)) +
## add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    ## custom bandwidth
    adjust = .5,
    ## adjust height
    width = .6,
    .width = 0,
    ## move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    ## remove outliers
    outlier.shape = NA ## `outlier.shape = NA` works as well
  ) +
  ## add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    ## draw jitter on the left
    side = "l",
    ## control range of jitter
    range_scale = .4,
    ## add some transparency
    alpha = .5,
    transformation = ggplot2::position_jitter(height = 0)

  ) +
  labs(y = "Decomposition constant (k)", x = "P treatment") +
  geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = n.treat, label = n.labels), nudge_x = 0, size = 6) +
     theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1))  +
    scale_x_discrete(labels=c("No P" = "No P", "P" = "+P"))

rn_n_wood

Code

ggsave("./output/nitrogen_decomposition_k.tiff", dpi = 300)

Download image

Code

fit_k.wood.n.no.n <- brm(sil.k.wood ~ n.treat + (1 | quadrat), data = k_vals,
    chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/k_wood_n_no_n_model",
    file_refit = "on_change")

Code

extended_summary(read.file = "./data/processed/k_wood_n_no_n_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

4.1 k_wood_n_no_n_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.wood ~ n.treat + (1 | quadrat)	20000	4	1	10000	7	0	30908.88	27135.72	1683655054

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
n.treatN	0.008	-0.414	0.426	1	30908.88	27135.72

5 Nutrient content by remaining litter mass

5.1 Nitrogen

Code

nutr$prom.hoja.reman.sc <- scale(nutr$prom.hoja.reman)

ggplot(data = nutr, aes(x = prom.hoja.reman, y = perc.n, color = treat)) +
    geom_point() + labs(x = "Remaining litter mass (% initial)", y = "Litter N concentration (%)",
    color = "Treatment") + scale_color_viridis_d(alpha = 0.5) + geom_smooth(method = "lm",
    se = FALSE) + scale_x_reverse()

Code

ggsave("./output/nitrogen_by_litter_mass.tiff", dpi = 300)

Download image

Code

fit_pern_by_rem_leave <- brm(perc.n ~ prom.hoja.reman.sc * treat +
    (1 | plot), data = nutr, chains = chains, family = gaussian(),
    iter = iters, control = list(adapt_delta = 0.99, max_treedepth = 15),
    cores = chains, prior = prior, file = "./data/processed/n_per_by_leave_rem_model",
    file_refit = "on_change")

Code

extended_summary(read.file = "./data/processed/n_per_by_leave_rem_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

5.2 n_per_by_leave_rem_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	perc.n ~ prom.hoja.reman.sc * treat + (1 | plot)	20000	4	1	10000	0	0	22712.6	26086.89	202030287

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
prom.hoja.reman.sc	-0.152	-0.265	-0.038	1	22712.60	27352.64
treatN	-0.161	-0.345	0.025	1	29404.33	26086.89
treatNP	-0.044	-0.224	0.138	1	29226.14	26624.22
treatP	-0.037	-0.218	0.145	1	28414.57	26316.85
prom.hoja.reman.sc:treatN	0.060	-0.095	0.212	1	27004.70	29358.59
prom.hoja.reman.sc:treatNP	-0.103	-0.269	0.064	1	28381.31	29425.85
prom.hoja.reman.sc:treatP	-0.013	-0.174	0.150	1	27640.20	30237.33

5.3 Phosphorus

Code

ggplot(data = sub.nutr, aes(x = prom.hoja.reman, y = perc.p, color = treat)) +
    geom_point() + labs(x = "Remaining litter mass (% initial)", y = "Litter P concentration (%)",
    color = "Treatment") + scale_color_viridis_d(alpha = 0.5) + geom_smooth(method = "lm",
    se = FALSE) + scale_x_reverse()

Code

ggsave("./output/phosphorus_by_litter_mass.tiff", dpi = 300)

Download image

Code

sub.nutr$prom.hoja.reman.sc <- scale(sub.nutr$prom.hoja.reman)

fit_perp_by_rem_leave <- brm(perc.p ~ prom.hoja.reman.sc * treat +
    (1 | plot), data = sub.nutr, chains = chains, family = gaussian(),
    iter = iters, control = list(adapt_delta = 0.99, max_treedepth = 15),
    cores = chains, prior = prior, file = "./data/processed/p_per_by_leave_rem_model",
    file_refit = "on_change")

Code

extended_summary(read.file = "./data/processed/p_per_by_leave_rem_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

5.4 p_per_by_leave_rem_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	perc.p ~ prom.hoja.reman.sc * treat + (1 | plot)	20000	4	1	10000	0	0	18469.62	19550.38	8818990

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
prom.hoja.reman.sc	0.001	-0.002	0.004	1	21082.94	25940.57
treatN	-0.005	-0.012	0.002	1	18785.39	19550.38
treatNP	0.000	-0.007	0.007	1	18469.62	20718.11
treatP	0.003	-0.004	0.011	1	18490.36	20587.53
prom.hoja.reman.sc:treatN	0.000	-0.004	0.005	1	24965.09	28102.51
prom.hoja.reman.sc:treatNP	-0.002	-0.006	0.003	1	25491.25	27830.25
prom.hoja.reman.sc:treatP	-0.001	-0.005	0.004	1	26931.48	27838.27

Takeaways

• Nitrogen percentage, but no Phosphorus percentage, increases along with remaining litter mass

 

6 Nutrient ratios

6.1 C:N

Code

# excluding plot 9
agg_cn <- aggregate(cn.mol.kg ~ colecta + treat + days, nutr, mean)

agg_cn$sd <- aggregate(cn.mol.kg ~ colecta + treat + days, nutr, sd)$cn.mol.kg

agg_cn$se <- aggregate(cn.mol.kg ~ colecta + treat + days, nutr, se)$cn.mol.kg


agg_cn$treat <- factor(agg_cn$treat, levels = c("C", "N", "P", "NP"))

pd <- position_dodge(15)

ggplot(agg_cn, aes(x = days, y = cn.mol.kg, color = treat)) + geom_point(size = 2,
    position = pd) + geom_errorbar(aes(ymax = cn.mol.kg + se, ymin = cn.mol.kg -
    se), width = 0, position = pd) + geom_line(size = 1.2, position = pd) +
    scale_color_viridis_d(alpha = 0.5) + labs(x = "Time (days)", y = "Litter C:N ratio",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_cn$days),
    labels = unique(agg_cn$days)) + theme(legend.position = c(0.9,
    0.8))

Code

ggsave("./output/cn_ratio_through_time.tiff", dpi = 300)

Download image

Code

fit.cn <- brm(cn.mol.kg ~ treat * days.sc + (1 | plot.f), data = nutr,
    chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/cn_model",
    file_refit = "on_change")

Code

extended_summary(read.file = "./data/processed/cn_model.rds", gsub.pattern = "b_treatment|b_",
    gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)

6.2 cn_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	cn.mol.kg ~ treat * days.sc + (1 | plot.f)	20000	4	1	10000	0	0	23616.14	27188.36	206225990

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
treatN	0.309	-1.373	2.015	1	33345.57	27188.36
treatNP	-0.189	-1.881	1.509	1	31802.40	27824.88
treatP	-0.182	-1.863	1.514	1	30821.46	28785.15
days.sc	-3.101	-4.204	-1.996	1	23616.14	27323.41
treatN:days.sc	0.751	-0.752	2.245	1	29469.69	30884.24
treatNP:days.sc	0.439	-1.062	1.942	1	29505.97	31955.60
treatP:days.sc	0.111	-1.391	1.614	1	29297.05	31633.29

6.3 N:P

Code

# excluding plot 9
agg_np <- aggregate(np.mol.kg ~ colecta + treat + days, sub.nutr,
    mean)

agg_np$sd <- aggregate(np.mol.kg ~ colecta + treat + days, sub.nutr,
    sd)$np.mol.kg

agg_np$se <- aggregate(np.mol.kg ~ colecta + treat + days, sub.nutr,
    se)$np.mol.kg


agg_np$treat <- factor(agg_np$treat, levels = c("C", "N", "P", "NP"))


pd <- position_dodge(15)

ggplot(agg_np, aes(x = days, y = np.mol.kg, color = treat)) + geom_point(size = 2,
    position = pd) + geom_errorbar(aes(ymax = np.mol.kg + se, ymin = np.mol.kg -
    se), width = 0, position = pd) + geom_line(size = 1.2, position = pd) +
    scale_color_viridis_d(alpha = 0.5) + labs(x = "Time (days)", y = "Litter N:P ratio",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_np$days),
    labels = unique(agg_np$days)) + theme(legend.position = c(0.2,
    0.8))

Code

ggsave("./output/np_ratio_through_time.tiff", dpi = 300)

Download image

Code

fit.np <- brm(np.mol.kg ~ treat * days.sc + (1 | plot.f), data = sub.nutr,
    chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/np_model",
    file_refit = "on_change")

Code

extended_summary(read.file = "./data/processed/np_model.rds", gsub.pattern = "b_treatment|b_",
    gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)

6.4 np_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 5.3)	np.mol.kg ~ treat * days.sc + (1 | plot.f)	20000	4	1	10000	0	0	23138.62	23571.24	226626333

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
treatN	1.211	-3.076	5.572	1	27662.79	23728.89
treatNP	-2.401	-6.601	1.896	1	26490.49	23571.24
treatP	-3.135	-7.696	1.500	1	27059.33	23906.00
days.sc	4.407	1.755	7.072	1	23138.62	26390.49
treatN:days.sc	-1.765	-5.391	1.885	1	26580.63	27306.41
treatNP:days.sc	-1.780	-5.403	1.901	1	27640.74	28593.42
treatP:days.sc	-1.415	-5.318	2.486	1	28596.79	29590.75

6.5 C:P

Code

# excluding plot 9
agg_cp <- aggregate(cp.mol.kg ~ colecta + treat + days, sub.nutr,
    mean)

agg_cp$sd <- aggregate(cp.mol.kg ~ colecta + treat + days, sub.nutr,
    sd)$cp.mol.kg

agg_cp$se <- aggregate(cp.mol.kg ~ colecta + treat + days, sub.nutr,
    se)$cp.mol.kg

agg_cp$treat <- factor(agg_cp$treat, levels = c("C", "N", "P", "NP"))


pd <- position_dodge(15)

ggplot(agg_cp, aes(x = days, y = cp.mol.kg, color = treat)) + geom_point(size = 2,
    position = pd) + geom_errorbar(aes(ymax = cp.mol.kg + se, ymin = cp.mol.kg -
    se), width = 0, position = pd) + geom_line(size = 1.2, position = pd) +
    scale_color_viridis_d(alpha = 0.5) + labs(x = "Time (days)", y = "Litter C:P ratio",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_cp$days),
    labels = unique(agg_cp$days)) + theme(legend.position = c(0.9,
    0.8))

Code

ggsave("./output/cp_ratio_through_time.tiff", dpi = 300)

Download image

Code

fit.cp <- brm(cp.mol.kg ~ treat * days.sc + (1 | plot.f), data = sub.nutr,
    chains = chains, family = gaussian(), iter = iters, cores = chains,
    prior = prior, file = "./data/processed/cp_model", file_refit = "on_change")

Code

extended_summary(read.file = "./data/processed/cp_model.rds", gsub.pattern = "b_treatment|b_",
    gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)

6.6 cp_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 106.3)	cp.mol.kg ~ treat * days.sc + (1 | plot.f)	20000	4	1	10000	0	0	30235.23	26541.32	1179254390

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
treatN	0.386	-19.206	20.050	1	30235.23	27580.58
treatNP	-0.159	-19.492	19.226	1	31461.97	26541.32
treatP	-0.300	-19.898	19.294	1	32550.30	27705.93
days.sc	-6.177	-21.408	9.060	1	30959.55	27830.18
treatN:days.sc	0.319	-17.923	18.305	1	32235.63	27315.59
treatNP:days.sc	-4.800	-23.038	13.518	1	32491.62	27863.26
treatP:days.sc	-2.648	-20.973	15.782	1	32066.82	27837.57

---

7 Combined plots

7.1 Remaining mass

Code

agg_rem_w$substrate <- "Wood"
agg_rem$substrate <- "Litter"

names(agg_rem) <- names(agg_rem_w) <- c("colecta", "trat", "days",
    "perc.rem", "sd", "se", "substrate")

agg_rem_pooled <- rbind(agg_rem, agg_rem_w)

ggplot(agg_rem_pooled, aes(x = days, y = perc.rem, color = trat)) +
    geom_point(size = 2, position = pd) + geom_errorbar(aes(ymax = perc.rem +
    se, ymin = perc.rem - se), width = 0, position = pd) + geom_line(size = 1.2,
    position = pd) + scale_color_viridis_d(alpha = 0.5, labels = c("Control",
    "+N", "+P", "+NP")) + labs(x = "Time (days)", y = "Remaining mass (%)",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_rem_pooled$days),
    labels = unique(agg_rem_pooled$days)) + facet_wrap(~substrate,
    nrow = 2)

Code

ggsave("./output/remaining_mass_combined.tiff", dpi = 300)

Download image

7.2 Litter content

Code

agg_p$nutrient <- "P"
agg_n$nutrient <- "N"

names(agg_n) <- names(agg_p) <- c("colecta", "treat", "days", "perc",
    "sd", "se", "nutrient")

agg_np2 <- rbind(agg_n, agg_p)

ggplot(agg_np2, aes(x = days, y = perc, color = treat)) + geom_point(size = 2,
    position = pd) + geom_errorbar(aes(ymax = perc + se, ymin = perc -
    se), width = 0, position = pd) + geom_line(size = 1.2, position = pd) +
    scale_color_viridis_d(alpha = 0.5, labels = c("Control", "+N",
        "+P", "+NP")) + labs(x = "Time (days)", y = "Litter nutrient content (% initial)",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_p$days),
    labels = unique(agg_p$days)) + facet_wrap(~nutrient, nrow = 2)

Code

ggsave("./output/litter_content_combined.tiff", dpi = 300)

Download image

7.3 Concentration

Code

agg_conc_p$nutrient <- "P"
agg_conc_n$nutrient <- "N"

names(agg_conc_n) <- names(agg_conc_p) <- c("colecta", "treat", "days",
    "perc", "sd", "se", "nutrient")

agg_conc_np <- rbind(agg_conc_n, agg_conc_p)

ggplot(agg_conc_np, aes(x = days, y = perc, color = treat)) + geom_point(size = 2,
    position = pd) + geom_errorbar(aes(ymax = perc + se, ymin = perc -
    se), width = 0, position = pd) + geom_line(size = 1.2, position = pd) +
    scale_color_viridis_d(alpha = 0.5, labels = c("Control", "+N",
        "+P", "+NP")) + labs(x = "Time (days)", y = "Litter nutrient concentration (%)",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_conc_p$days),
    labels = unique(agg_conc_p$days)) + facet_wrap(~nutrient, nrow = 2,
    scales = "free_y")

Code

ggsave("./output/concentration_combined.tiff", dpi = 300)

Download image

7.4 Nutrient ratios

Code

agg_cn$nutrients <- "C:N"
agg_cp$nutrients <- "C:P"
agg_np$nutrients <- "N:P"

names(agg_cn) <- names(agg_cp) <- names(agg_np) <- c("colecta", "treat",
    "days", "mol.kg", "sd", "se", "nutrients")

agg_ratios <- rbind(agg_cn[, names(agg_cn)], agg_cp[, names(agg_cn)],
    agg_np[, names(agg_cn)])

ggplot(agg_ratios, aes(x = days, y = mol.kg, color = treat)) + geom_point(size = 2,
    position = pd) + geom_errorbar(aes(ymax = mol.kg + se, ymin = mol.kg -
    se), width = 0, position = pd) + geom_line(size = 1.2, position = pd) +
    scale_color_viridis_d(alpha = 0.5, labels = c("Control", "+N",
        "+P", "+NP")) + labs(x = "Time (days)", y = "Litter nutrient ratios",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_ratios$days),
    labels = unique(agg_ratios$days)) + facet_wrap(~nutrients, nrow = 3,
    scales = "free_y")

Code

ggsave("./output/nutrient_ratios_combined.tiff", dpi = 300, width = 10,
    height = 12)

Download image

7.5 Decomposition constant

Code

rn_n_wood_dat <- rn_n_wood$data[, c("n.treat", "sil.k.wood")]
rn_n_litter_dat <- rn_n_litter$data[, c("n.treat", "sil.k.litter")]
rn_p_wood_dat <- rn_p_wood$data[, c("p.treat", "sil.k.wood")]
rn_p_litter_dat <- rn_p_litter$data[, c("p.treat", "sil.k.litter")]

names(rn_p_wood_dat) <- names(rn_n_wood_dat)  <- names(rn_p_litter_dat) <- names(rn_n_litter_dat) <- c("treatment", "k")
rn_p_wood_dat$substrate <- rn_n_wood_dat$substrate <- "wood"
rn_p_litter_dat$substrate <- rn_n_litter_dat$substrate <- "litter"

rn_n_wood_dat$nutrient <- rn_n_litter_dat$nutrient <- "N"
rn_p_litter_dat$nutrient <- rn_p_wood_dat$nutrient <- "P"

rn_k_dat <- rbind(rn_n_wood_dat, rn_p_wood_dat, rn_n_litter_dat, rn_p_litter_dat)                             
                             
# composed box plot
ggplot(rn_k_dat, aes(x = treatment, y = k)) +
## add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    ## custom bandwidth
    adjust = .5,
    ## adjust height
    width = .6,
    .width = 0,
    ## move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    ## remove outliers
    outlier.shape = NA ## `outlier.shape = NA` works as well
  ) +
  ## add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    ## draw jitter on the left
    side = "l",
    ## control range of jitter
    range_scale = .4,
    ## add some transparency
    alpha = .5,
    transformation = ggplot2::position_jitter(height = 0)

  ) +
  labs(y = "Decomposition constant (k)", x = "P treatment") +
  # geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = n.treat, label = n.labels), nudge_x = 0, size = 6) +
     theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1))  +
    scale_x_discrete(labels=c("No P" = "No P", "P" = "+P")) +
    facet_grid(~ substrate)

Code

ggsave("./output/decomposition_k_combined_v1.tiff", dpi = 300, width = 10, height = 12)

Download image

Code

# composed box plot
ggplot(rn_k_dat, aes(x = treatment, y = k)) +
## add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    ## custom bandwidth
    adjust = .5,
    ## adjust height
    width = .6,
    .width = 0,
    ## move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    ## remove outliers
    outlier.shape = NA ## `outlier.shape = NA` works as well
  ) +
  ## add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    ## draw jitter on the left
    side = "l",
    ## control range of jitter
    range_scale = .4,
    ## add some transparency
    alpha = .5,
    transformation = ggplot2::position_jitter(height = 0)

  ) +
  labs(y = "Decomposition constant (k)", x = "P treatment") +
  # geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = n.treat, label = n.labels), nudge_x = 0, size = 6) +
     theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1))  +
    scale_x_discrete(labels=c("No P" = "No P", "P" = "+P")) +
    facet_grid(substrate ~ nutrient, scales = "free_x")

Code

ggsave("./output/decomposition_k_combined_v12.tiff", dpi = 300, width = 10, height = 12)

Download image

8 Combined model diagnostics

Code

check_rds_fits(path = "./data/processed", html = TRUE)

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	cn.mol.kg ~ treat * days.sc + (1 | plot.f)	20000	4	1	10000	0	0	11186.547	18465.237	206225990
2	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 106.3)	cp.mol.kg ~ treat * days.sc + (1 | plot.f)	20000	4	1	10000	0	0	3800.096	3310.208	1179254390
3	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.litter ~ Treatment + (1 | quadrat)	20000	4	1	10000	5	0	10254.769	11806.319	224004760
4	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.litter ~ n.treat + (1 | quadrat)	20000	4	1	10000	17	0	8997.572	8327.074	956380015
5	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.litter ~ p.treat + (1 | quadrat)	20000	4	1	10000	8	0	10806.876	12607.636	1236972411
6	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.wood ~ n.treat + (1 | quadrat)	20000	4	1	10000	7	0	8165.879	9671.253	1683655054
7	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.wood ~ p.treat + (1 | quadrat)	20000	4	1	10000	4	0	8190.692	11336.204	326631910
8	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.wood ~ Treatment + (1 | quadrat)	80000	4	1	40000	23	0	37124.811	54565.753	528819192
9	b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20)	litter.n.content.prop.initial ~ treat * days.sc + (1 | plot.f)	20000	4	1	10000	0	0	12453.012	17014.307	272268972
10	b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20)	prop.n ~ treat * days.sc + (1 | plot.f)	20000	4	1	10000	0	0	14603.282	16880.616	757288462
11	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	litter.p.content.prop.initial ~ treat * days.sc + (1 | plot.f)	20000	4	1	10000	0	0	8457.580	15360.879	2063397055
12	b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20)	prop.p ~ treat * days.sc + (1 | plot.f)	20000	4	1	10000	0	0	16204.903	17372.163	1797840914
13	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	perc.n ~ prom.hoja.reman.sc * treat + (1 | plot)	20000	4	1	10000	0	0	12014.609	18513.789	202030287
14	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 5.3)	np.mol.kg ~ treat * days.sc + (1 | plot.f)	20000	4	1	10000	0	0	11291.391	17616.689	226626333
15	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	perc.p ~ prom.hoja.reman.sc * treat + (1 | plot)	20000	4	1	10000	0	0	7081.653	10836.794	8818990
16	b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20)	prop.litter.rem ~ trat * days.sc + (1 | plot.f)	20000	4	1	10000	0	0	13040.090	18549.848	1631599382
17	b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20)	prop.wood.rem ~ trat * days.sc + (1 | plot.f)	20000	4	1	10000	0	0	11933.976	17029.142	1828689869

 

Session information

R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3

locale:
 [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=pt_BR.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=pt_BR.UTF-8   
 [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ggdist_3.2.0      brms_2.18.0       Rcpp_1.0.9        viridis_0.6.2    
 [5] viridisLite_0.4.1 ggridges_0.5.4    posterior_1.3.1   ggrepel_0.9.1    
 [9] cowplot_1.1.1     tidybayes_3.0.1   ggplot2_3.4.0     readxl_1.3.1     
[13] knitr_1.42        kableExtra_1.3.4 

loaded via a namespace (and not attached):
  [1] backports_1.4.1      systemfonts_1.0.4    plyr_1.8.7          
  [4] igraph_1.3.5         splines_4.1.0        svUnit_1.0.6        
  [7] crosstalk_1.2.0      TH.data_1.1-0        rstantools_2.2.0    
 [10] inline_0.3.19        digest_0.6.31        htmltools_0.5.4     
 [13] fansi_1.0.3          magrittr_2.0.3       checkmate_2.1.0     
 [16] remotes_2.4.2        RcppParallel_5.1.5   matrixStats_0.62.0  
 [19] xts_0.12.2           sandwich_3.0-1       svglite_2.1.0       
 [22] prettyunits_1.1.1    colorspace_2.0-3     rvest_1.0.3         
 [25] textshaping_0.3.5    xfun_0.36            dplyr_1.0.10        
 [28] callr_3.7.3          crayon_1.5.2         jsonlite_1.8.4      
 [31] lme4_1.1-27.1        survival_3.2-11      zoo_1.8-11          
 [34] glue_1.6.2           gtable_0.3.1         emmeans_1.8.1-1     
 [37] webshot_0.5.4        distributional_0.3.1 pkgbuild_1.4.0      
 [40] rstan_2.21.7         abind_1.4-5          scales_1.2.1        
 [43] mvtnorm_1.1-3        DBI_1.1.1            xaringanExtra_0.7.0 
 [46] miniUI_0.1.1.1       xtable_1.8-4         stats4_4.1.0        
 [49] StanHeaders_2.21.0-7 DT_0.26              htmlwidgets_1.5.4   
 [52] httr_1.4.4           threejs_0.3.3        arrayhelpers_1.1-0  
 [55] ellipsis_0.3.2       pkgconfig_2.0.3      loo_2.4.1.9000      
 [58] farver_2.1.1         utf8_1.2.2           labeling_0.4.2      
 [61] tidyselect_1.2.0     rlang_1.0.6          reshape2_1.4.4      
 [64] later_1.3.0          munsell_0.5.0        cellranger_1.1.0    
 [67] tools_4.1.0          cli_3.6.0            generics_0.1.3      
 [70] evaluate_0.20        stringr_1.5.0        fastmap_1.1.0       
 [73] ragg_1.1.3           yaml_2.3.7           processx_3.8.0      
 [76] purrr_1.0.0          packrat_0.9.0        pbapply_1.6-0       
 [79] nlme_3.1-152         mime_0.12            projpred_2.0.2      
 [82] formatR_1.11         xml2_1.3.3           compiler_4.1.0      
 [85] bayesplot_1.9.0      shinythemes_1.2.0    rstudioapi_0.14     
 [88] gamm4_0.2-6          sketchy_1.0.2        tibble_3.1.8        
 [91] stringi_1.7.12       ps_1.7.2             Brobdingnag_1.2-9   
 [94] lattice_0.20-44      Matrix_1.5-1         nloptr_1.2.2.2      
 [97] markdown_1.3         shinyjs_2.1.0        tensorA_0.36.2      
[100] vctrs_0.5.2          pillar_1.8.1         lifecycle_1.0.3     
[103] bridgesampling_1.1-2 estimability_1.4.1   httpuv_1.6.6        
[106] R6_2.5.1             promises_1.2.0.1     gridExtra_2.3       
[109] codetools_0.2-18     boot_1.3-28          colourpicker_1.2.0  
[112] MASS_7.3-54          gtools_3.9.3         assertthat_0.2.1    
[115] rprojroot_2.0.3      withr_2.5.0          shinystan_2.6.0     
[118] multcomp_1.4-17      mgcv_1.8-36          parallel_4.1.0      
[121] gghalves_0.1.3       grid_4.1.0           minqa_1.2.4         
[124] tidyr_1.1.3          coda_0.19-4          rmarkdown_2.20      
[127] shiny_1.7.3          base64enc_0.1-3      dygraphs_1.1.1.6